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INTRODUCTION 

Most  cartridge-based  munitions  use  some  form  of  granular  solid  propellant  (fig.  1),  in  which 
the  deflagration  process  is  modeled  as  a  regression  of  grain  surfaces  along  their  respective  surface 
normal  (Piobert’s  Law  of  Burning,  1839).  The  propellant  gas  generation  rate  is  well  known  to  be 
dependent  on  both  local  pressure  and  the  granular  surface  area.  Control  over  the  shape  of  solid 
propellant  grain  allows  for  a  high  degree  of  control  over  the  propellant  burn  rate,  burn  time,  and 
amount  of  generated  propellant  gases  that  directly  impact  the  thrust  versus  time  profile  of  the  given 
system  (ref.  1).  Generally,  grains  with  high  progressivity,  or  grains  whose  burning  surface  area 
increase  as  the  grain  is  consumed,  are  desired  so  that  the  generation  rate  of  propellant  gases 
increases  during  the  munition  launch  cycle.  Some  typical  approaches  to  increase  propellant  grain 
surface  area  are  the  addition  of  slotted/perforated  features  to  the  grain  geometry  (fig.  2). 


Figure  1 

Cartridge-based  munition  containing  single-perforation  cylindrical  propellant  grains 


(a)  (b) 

M549  aft  rocket  motor  XM 1 1 28  rocket  motor 

Figure  2 

Rocket  motor  geometries  with  slotted/perforated  features  to  increase  burn  surface  area 
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A  common  mathematical  model  that  relates  grain  linear  surface  regression  rate  dx/dt  to 
current  local  pressure  P;(t)  is  Vieille’s  Law  of  1893,  shown  in  equation  1 

(1) 

where  x  equals  x(t)  is  the  linear  surface  regression  at  time  t,  and  A  and  n  are  constants  that  depend 
on  the  propellant  composition.  For  simple  propellant  grain  geometries,  such  as  the  single¬ 
perforation  [perforation  (cylinder  with  concentric  hole)]  grain,  mathematical  expressions  for  the 
unburnt  grain  mass  mg(t)  and  gas  generation  rate  dmg/dt  at  a  given  time  t  can  be  readily 
determined  from  the  initial  grain  geometry  as  shown  in  figure  3.  In  all  forthcoming  sections,  the 
instantaneous  linear  surface  regression  x(t)  has  been  written  as  x  for  convenience. 


TT  . 


mg(t)  =  Pg*Vg(t)  =  pg'-[(D  -  2  *x)2  -  (P  +  2  *x)2](L  -2  *  x) 


dmg(t) 


dx 

drrig^t)  dmg(t)  dx 


,n  i  l  j  TC 

9  =  -pg-(D  +  P)(D  -  P  +  2L  —  8  *  x) 


dt 


dx  dt 


-AP?(t)pg-(D  +  P)(D  -  P  +  2L  -  8  *  x) 


Figure  3 

Surface  regression  diagram  and  burn  equations  for  single-perforation  grain  with  initial  outer  diameter 

D,  inner  diameter  P,  and  length  L  [not  shown  (ref.  2)] 

The  advent  of  additive  manufacturing  brings  forth  a  new  advantage  in  solid  propellant 
technology  -  the  ability  to  create  grains  of  arbitrary  shape  with  complex  internal  structures  that  allow 
for  fine-tuning  of  the  burning  surface  area  over  time,  and  thus,  greater  control  of  propellant  burn  rate. 
A  consequence  of  these  complex,  additive  geometries,  however,  is  that  analytical  expressions  and 
prediction  of  burn  rates  from  surface  regression  can  become  infeasible  or  impossible  to  obtain.  In 
order  to  facilitate  additive  propellant  grain  design  optimization  without  requiring  physical  printing  and 
costly  burn  testing,  a  numerical  tool  predicting  temporal  burn  regression/rate  curves  from  initial 
computer  aided  design  (CAD)  geometries  is  required. 

To  this  end,  the  U.S.  Army  Armament  Research,  Development  and  Engineering  Center 
(ARDEC),  Picatinny  Arsenal,  NJ,  developed  the  ARDEC  grain  evaluation  software  (AGES),  a 
numerical  tool  to  calculate  burn  regression  and  burn  rate  curves  for  arbitrary  solid  propellant 
geometry  through  use  of  phase-field  interface  tracking  of  the  burning  grain  surface.  The  AGES  uses 
a  Eulerian  volume  of  fluid  approach  to  “track”  the  exterior  surface  of  the  grain  over  time.  The 
numerical  approach  and  development  of  AGES  is  outlined  in  this  report,  and  numerical  results  are 
compared  to  analytical  expressions  for  common  propellant  grain  geometries  as  validation.  Finally, 
AGES  is  applied  to  prototypical  additively-manufactured  grain  geometry  with  complex  internal 
geometries  to  showcase  the  effectiveness  of  this  numerical  approach. 
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PHASE-FIELD  EQUATION 

Numerical  tracking  of  interfaces  using  an  Eulerian  (fixed-grid)  method  is  commonly 
implemented  in  situations  where  explicit  tracking  of  exact  interface  locations/coordinates  are 
mathematically  difficult  (ref.  3).  Examples  of  applications  involve  multiphase  flows,  co-extrusion,  and 
phase-change  problems.  The  general  phase-field  method  for  sharp  interface  tracking  is 
implemented  in  AGES.  Using  this  method,  moving  interfaces  are  tracked  via  the  solution  of  the 
general  phase-field  interface  advection  equation 

^  +  u-V0  =  O  (2) 

at 

where  0  is  the  phase-field  and  u  is  the  interfacial  velocity.  Various  kernel  functions  can  be  used  to 
describe  the  variation  of  the  phase-field  0  normal  to  the  interface.  A  common  kernel  function  for  0  is 
given  by  equation  3  (ref.  4) 


where  x'  is  the  ratio  between  the  normal  distance  of  a  point  in  the  phase-field  from  the  interface  n 
and  the  width  of  the  hyperbolic  tangent  profile  W.  The  variation  of  0  with  x'  is  shown  in  figure  4  for 
the  hyperbolic  tangent  kernel.  As  the  distance  from  the  interface  reaches  a  critical  value,  the  phase- 
field  value  approaches  a  constant  -  a  distinct  advantage  from  the  level-set  method  where  the  signed 
distance  always  increases.  At  the  interface,  the  value  of  0  is  exactly  0,  while  at  locations  sufficiently 
far  away  from  the  interface,  |0|  equals  1. 


Figure  4 

Variation  of  phase-field  value  0  around  interface  (*’  equals  0) 


For  an  interface  that  has  no  curvature-driven  motion,  equation  2  becomes 

|  V0 1 V  ■ 


d0 

+  fl|V0|  =  b 


^.0(i zTl 


Vz0  +  - 


w2 


MV0I/ 


(4) 


where  a  is  the  interface  velocity  and  b  is  a  purely  numerical  parameter  that  controls  the  relaxation 
and  smoothing  of  the  interface.  The  terms  on  the  left-hand  side  represent  temporal  and  spatial 
changes  in  phase-field  values,  respectively.  The  first  two  terms  on  the  right-hand  side  represent 
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interface  curvature,  and  the  third  term  is  a  curvature  counter  term  that  mathematically  cancels  the 
interface  curvature  terms  at  the  leading  order,  while  maintaining  the  hyperbolic  tangent  profile  at  the 
interface  (ref.  5).  For  a  more  detailed  derivation  and  explanation  of  the  governing  equation,  see 
reference  4. 


NUMERICAL  DISCRETIZATION 

Numerical  implementation  of  equation  4  can  be  achieved  using  finite-differencing.  Equation 
4  is  discretized  on  a  three-dimensional  rectangular  grid  with  equal  spacing  in  all  three  ordinate 
directions  (i.e.,  Ax  equals  Ay  equals  A z).  This  is  critical  for  discretization  of  the  Laplacian  term  V20, 
which  can  be  discretized  using  the  highly-stable  19-point  stencil  (ref.  6) 

^20OOO  —  [2(0+00  +  0-00  +  00+0  +  00-0  +  000+  +  000-) 

+  (0++o  +  0-+o  +  0+-o  +  0 — o  +  0+o+  +  0+o-  +  0-o+  +  0-o-  +  0o++  +  0  o+-  (5) 

+  0o — v  +  0o — )  —  240ooo]/  (6Ax2) 

where  the  notation  0+o_  indicates  0i+1,;,fc-1  in  traditional  finite-differencing  notation.  Although  the 
choice  was  made  to  implement  a  higher-order  approximation  to  the  Laplacian,  the  standard  7-point 
stencil  for  the  Laplacian  can  be  implemented  without  significant  error  for  most  cases.  The  gradient 
V0  is  discretized,  using  central  differencing. 


Finally,  the  curvature  term  is  discretized  in  the  following  manner 


where  0Cx+  and  0Cx_  can  be  expressed  as 

4>Cx+ 


0+oo  —  0ooo 


and  the  remaining  terms  inside  the  square  bracket  can  be  expressed  in  a  similar  manner. 

The  temporal  derivative  in  equation  4  is  discretized  using  first-order  forward-differencing  with 
time  step  At.  The  numerical  parameter  W,  corresponding  to  the  width  of  the  hyperbolic  tangent 
profile  in  figure  4,  must  be  defined  small  enough  to  resolve  the  0  profile  numerically.  To  avoid 
overlapping  of  0  profiles  for  a  curved  interface,  W  should  be  kept  as  small  as  possible,  although  it  is 
absolutely  necessary  that  W  >  Ax.  Parametric  studies  on  W  have  shown  that  a  value  of 
W  equals  2Ax  is  suitable  for  most  situations  (ref.  4). 
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Finally,  the  numerical  parameter  b,  which  helps  control  problem  stability  and  convergence, 
should  adhere  to  the  following  relation 


3Ax2  ,, 

b<im^~2Cr)  {i 

as  described  in  reference  4,  where  the  Courant  number  Cr  equals  aAt/Ax.  Typically,  a  value  of 
Cr  equals  0.1  provides  sufficient  accuracy  without  great  increase  in  computational  cost.  The 
interface  velocity  a,  in  a  sense,  is  therefore  unimportant  except  to  determine  the  necessary  At  to 
achieve  the  desired  Cr  for  a  given  grid  spacing.  Therefore,  for  all  simulations  forthcoming,  the 
interface  velocity  was  set  at  a  constant  0.5,  and  thus,  the  relationship  between  At  and  Ax  becomes 
At  equals  0.05Ax.  Using  this  relation,  equation  9  reduces  to  b  <  4.8Ax. 


PHASE-FIELD  INITIALIZATION 

Initialization  of  the  numerical  phase-field  stems  from  the  volume  fraction  of  grain  material  that 
exists  in  each  cell  in  the  computational  domain  at  t  equals  0.  Consider  the  two-dimensional  (2D), 
five  by  five  computational  grid  in  figure  5a.  For  the  circular  grain  shown  (blue),  the  volume  fraction  of 
grain  material  existing  in  each  computational  cell  is  listed,  ranging  from  0.0  (no  grain  material  exists 
in  the  cell)  to  1 .0  (the  entire  cell  is  comprised  of  grain  material). 
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Figure  5 

Phase-field  initialization 

To  convert  cell  volume  fraction  to  phase-field  values,  the  following  relationship  is  applied 

(l>Cell  =  2*VFcell-l  (10) 

Using  this  relationship,  a  cell  that  is  fully  filled  with  grain  material  will  have  (pceU  equals  1.0, 
and  a  cell  that  is  fully  void  will  have  0ceii  equals  -  1.0.  These  values  echo  the  concept  illustrated  in 
figure  4,  as  cells  that  are  either  fully  filled  or  fully  empty  do  not  intersect  with  the  boundary  interface. 
Phase-field  values  corresponding  to  the  volume  fractions  of  figure  5a  are  shown  in  figure  5b. 
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The  linear  relationship  of  equation  10  is  further  justified  when  considering  the  special  case  of 
VFCeii  equals  0.50,  where  exactly  half  of  the  cell  is  filled  with  material.  In  the  limit  of  zero  boundary 
curvature  (i.e.,  when  Ax  becomes  sufficiently  small),  the  boundary  interface  of  the  grain  will  pass 
directly  through  the  cell  center.  As  shown  in  figure  4,  this  situation  corresponds  to  0ceU  equals  0.0, 
or  zero  signed  distance  between  the  cell  center  and  the  boundary  interface.  Substitution  of 
VFCeii  equals  0.50  into  equation  10  produces  the  necessary  (pceU  equals  0.0. 


U.S.  ARMY  ARMAMENT  RESEARCH,  DEVELOPMENT  AND  ENGINEERING  CENTER  GRAIN 
PREPROCESSING  WITH  ABAQUS/FINITE  ELEMENT  ANALYSIS  (FEA) 

In  order  to  generate  the  phase-field  inputs  to  AGES,  the  built-in  volume  fraction  tool  in 
ABAQUS/FEA,  a  commercial  finite  element  package,  is  appropriated.  Solid  grain  geometry  can  be 
imported  into  ABAQUS  in  a  variety  of  different  file  formats.  In  conjunction  with  AGES,  multiple 
Python  scripts  were  developed  to  aid  in  the  creation  of  the  computational  domain,  which  is 
generated  using  Eulerian  elements.  Once  the  propellant  geometry  is  imported  into  ABAQUS,  the 
user  can  set  three  variables:  (1)  the  number  of  computational  elements  that  will  discretize  the 
shortest  side  of  the  computational  domain,  (2)  the  domain  extension  past  the  minimum  geometric 
bounds  of  the  propellant  grain,  and  (3)  the  accuracy  of  the  volume  fraction  tool  (low,  medium,  or 
high).  Figure  6  shows  a  typical  Eulerian  computational  domain  generated  in  ABAQUS  for  a  7- 
perforation  hex  propellant  grain.  Once  the  computational  Eulerian  domain  is  generated,  the  built-in 
volume  fraction  tool  can  be  run  to  get  the  fraction  of  each  Eulerian  cell  that  overlaps  the  solid 
geometry.  This  volume  fraction  is  then  converted  into  phase-field  values  using  equation  10  and 
exported  to  a  file  for  processing  in  the  ARDEC  grain  evaluation  software. 


(a) 

Eulerian  domain  with  computational 
mesh  visualized  as  generated 
in  ABAQUS 


(b) 

Eulerian  domain  with  propellant 
grain  geometry  fully  enclosed 


Figure  6 

Typical  Eulerian  computational  domain  generated  in  ABAQUS  for  a  7-perforation  hex  propellant 

grain 
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RESULTS  AND  DISCUSSION 


Validation,  Spherical  “Ball”  Grain 


As  an  initial  validation  case,  consider  a  spherical  “ball”  grain  of  an  initial  unit  diameter 
(Dg  equals  1.0).  For  spherical  geometry,  analytical  expressions  for  the  grain  volume  Vg  and 
instantaneous  change  in  grain  volume  dVg/dx  with  respect  to  surface  regression  x  can  be  written  as 


=  z(Ds-  2xf 


dVg(x) 

dx 


—  —n(pg  —  2x)2 


(11) 


The  quantity-  dVg/dx,  which  can  be  thought  of  as  an  instantaneous  change  in  volume 
consumption,  relates  directly  to  the  amount  of  propellant  gas  generated  during  the  burn  process. 

Figure  7a  compares  results  from  AGES  for  grain  volume  Vg  and  change  in  volume 
consumption  -dVg/dx  versus  burn  distance  x  for  the  unit  diameter  spherical  grain  to  the  exact 
solutions  of  equation  1 1 .  The  computational  domain  for  this  geometry  was  cubic  in  shape  with  a  side 
length  of  1 .02  (in  order  to  fully  encapsulate  the  unit  diameter  sphere).  The  domain  was  discretized 
using  either  N  equals  50,  100,  or  200  elements  per  side,  with  the  total  number  of  elements  in  the 
domain  equal  to  N3.  For  this  spherical  geometry,  total  grain  burnout  should  occur  at  x  equals  Dg/2 
equals  0.50. 


(a) 

Grain  volume  and  change  in  volume  consumption 


(b) 

Relative  percentage  error  in  grain  volume 
versus  burn  distance  for  spherical  ball  grain 


Figure  7 

Comparisons  and  results  generated  by  AGES 

Figure  7a  shows  that  numerical  results  generated  by  AGES  conform  well  to  the  analytical 
expressions  in  equation  11.  Errors  in  -dVg/dx  at  initial  small  burn  distances  are  witnessed  due  to 
the  relaxation  of  the  phase-field,  as  well  as  the  faceted  approximation  of  the  smooth  spherical  ball 
surface  with  cubic  elements.  As  seen  in  figure7b,  the  relative  percentage  error  between  the 
numerical  and  analytical  approximations  of  grain  volume  Vg  decrease  with  increasing  element  count. 
Errors  are  less  than  1  %  for  N  equals  200  until  x  equals  0.387,  at  which  point  only  ~1 .2%  of  the  initial 
grain  volume  remains  unburnt.  As  the  sphere  volume  becomes  significantly  small,  the  initial  element 
discretization  becomes  unsuitable  to  properly  capture  grain  regression,  leading  to  larger  errors.  For 
N  equals  200,  errors  of  >10%  occur  at  x  >  0.485,  at  which  point  the  percent  of  unburnt  grain  volume 
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is  less  than  0.1%.  Total  burnout  for  N  equals  200  occurs  at  x  equals  0.496,  which  is  within  1%  of  the 
true  value. 

As  previously  discussed,  the  volume  fraction  tool  in  ABAQUS/FEA  used  to  generate  the  initial 
phase-field  has  three  accuracy  methods:  (1)  low,  (2)  medium,  and  (3)  high  accuracy.  Figure  8 
shows  a  comparison  of  grain  volume  Vg  and  relative  error  in  Vg  versus  burn  distance  x  for  the 
spherical  grain  for  all  three  volume  fraction  tool  accuracies  using  grid  size  of  N  equals  200.  As 
expected,  high  accuracy  produces  the  lowest  overall  errors  as  burn  distance  increases,  although  for 
a  brief  period  after  initial  phase-field  relaxation,  medium  accuracy  produces  more  accurate  solutions. 
The  initial  grain  volume  calculated  using  the  three  accuracies  are  0.5178,  0.5196,  and  0.5226, 
respectively.  As  compared  to  the  true  initial  grain  volume  of  0.5236,  it  is  clear  that  running  the 
volume  fraction  tool  at  high  accuracy  is  preferred  to  minimize  numerical  error. 


Figure  8 

Grain  volume  and  relative  percentage  error  versus  burn  distance  for  spherical  ball  grain  generated 

with  low,  medium,  or  high  volume  fraction  accuracy 

Single-Perforation  Grain,  Superposition  Approach 

A  common  propellant  grain  shape  used  in  various  munitions  is  the  single-perforation,  which  is 
a  cylindrical  grain  of  diameter  Dg  and  length  Lg  with  a  concentric  hole  of  diameter  Pg  passing  through 
the  entire  length.  For  such  a  grain,  analytical  expressions  for  Vg  and  dVg/dx  for  a  given  burn 
distance  x  are  easily  calculated: 

\  ds  -  2 *)  [(»„  -  2 xf  -  (P„  +  2xf] 

=  [(D»  ”  2* )2  -  to  +  2*)2]  +  i  to  -  2*)H(»„  -  2x)  -  4 (Ps  +  2*)]  <12> 

Figure  9  compares  grain  volume  Vg  and  change  in  volume  consumption  -dVg/dx  versus  burn 
distance  x  for  a  single-perforation  propellant  grain  with  initial  Dg  equals  1.0,  Lg  equals  2.0,  and  Pg 
equals  0.5.  For  such  a  grain,  total  burnout  will  occur  at  x  equals  0.125.  The  computational  domain 
has  dimensions  x  equals  1 .02,  y  equals  1 .02,  and  z  equals  2.04  to  fully  enclose  the  grain. 
Discretization  of  the  domain  is  performed  with  three  element  densities: 

Nx  equals  Ny  equals  0.5 Nz  equals  101 , 201 ,  and  301 .  As  seen  for  the  spherical  grain  in  figure  7a, 
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errors  in  -dVg/dx  at  small  x  occur  due  to  initial  phase-field  relaxation  and  numerical  faceting.  After 
initial  relaxation,  numerical  results  conform  well  to  the  analytical  expressions  in  equation  12  until 
burn  distance  becomes  large.  When  10%  of  the  grain  volume  remains  at  x  equals  0.075,  the  relative 
percentage  errors  in  Vg  are  26.1%,  8.87%,  and  4.19%  for  the  three  element  densities,  respectively. 
Further  grid  refinement,  while  computationally  expensive,  would  improve  discrepancies  between 
analytical  and  numerical  results. 


Single  Perf  Grain  Regression 


1.20  12.00 


Burn  Distance,  x 


- Volume,  Nx  =  101  - Volume,  Nx  =  201  - Volume,  Nx  =  301  - Volume,  Analytical 

- -dV/dx.  Nx  =  101 - -dV/dx,  Nx  =  201 - -dV/dx,  Nx  =  301 - -dV/dx,  Analytical 


Figure  9 

Grain  volume  and  change  in  volume  consumption  versus  burn  distance  for  single-perforation  grain 

As  the  single-perforation  grain  approaches  burnout,  large  numerical  errors  occur,  which 
result  in  largely  inaccurate  burnout  distances  of  0.182,  0.162,  and  0.151  for  the  three  element 
densities,  respectively.  These  large  errors  stem  from  the  fact  that,  for  the  single-perforation  grain, 
there  are  two  interfaces  that  are  being  tracked:  (1)  the  outer  boundary  of  the  grain,  and  (2)  the 
interior  perforated  boundary.  When  these  two  interfaces  approach  one  another  during  burn, 
significant  diffusion  between  the  two  interfaces  occurs.  Additionally,  the  initial  grid  discretization 
becomes  too  small  to  accurately  capture  the  volume  in  the  small  “sliver”  of  remaining  grain 
geometry.  Numerical  errors  due  to  interface  diffusion  will  manifest  much  more  greatly  for  grains  with 
multiple  interacting  perforations,  such  as  the  common  7  or  19-perforation  rod,  and  thus,  a  method  to 
mitigate  these  types  of  errors  is  desired. 

One  method  to  mitigate  interface  diffusion  errors  is  by  applying  a  superposition  technique. 
Consider  the  2D  approximation  to  the  single-perforation  grain  in  figure  10.  Single-perforation  grain 
regression  can  be  effectively  described  as  the  superposition  of  two  solid  cylindrical  grains: 
regression  of  a  solid  grain  with  diameter  Dg  minus  the  progression  of  a  solid  grain  of  diameter  Pg.  By 
treating  a  single-perforation  grain  as  the  superposition  subtraction  of  two  solid  grains,  the  issue  of 
multiple  interface  interaction/diffusion  is  entirely  eliminated. 
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Illustration  of  superposition  for  2D  single-perforation  grain 

Using  the  superposition  in  figure  10,  two  phase-fields  are  initially  generated:  (1)  <pe,  which 
captures  the  external  boundary  defined  by  Dg,  and  (2)  0 £,  which  captures  the  internal  boundary 
defined  by  Pg.  The  internal  phase-field  0*  is  initialized  in  the  opposite  manner  of  0e  (i.e., 
cf>ceii  equals  -  1  inside  the  boundary,  and  (pceU  equals  1  outside  the  boundary)  so  that  the  internal 
boundary  will  grow  using  the  phase-field  equation  formulation.  The  volume  of  solid  propellant  in  a 
single  computational  cell,  therefore,  can  be  calculated  as: 

Vceli=AxAyAz[VFe-VFi]  (13) 


where 


VFe=\[4>e  +  1] 
VFi=l[l-4>i] 


(14) 


This  methodology  can  be  extended  to  grains  with  internal  features  or  grains  with  multiple 
perforations  that  may  interact  with  each  other  during  grain  regression. 

In  AGES,  a  scan  over  the  elements  in  the  computational  domain  is  performed  to  determine 
the  number  of  distinct  internal  features/cavities  that  exist  in  the  grain  geometry.  These  features  can 
either  be  set  to  all  burn  immediately,  or  to  burn  only  when  a  currently  burning  grain/feature  boundary 
interacts  with  them  for  the  first  time.  It  should  be  mentioned  that  for  each  internal  feature  or  cavity  in 
the  grain,  the  necessary  computational  memory  resources  will  approximately  double  since  a  full 
phase-field  for  the  full  computational  domain  needs  to  be  stored  for  each  progressing/regressing 
feature. 


For  grains  with  fully  enclosed  internal  features,  it  is  straightforward  to  determine  which 
computational  cells  encompass  those  features  through  a  simple  scanning  algorithm.  However,  for 
grains  with  perforations  through  the  entire  grain  (such  as  the  single-perforation),  it  becomes  more 
difficult  to  separate  the  features  numerically.  One  method  of  ensuring  that  the  perforation  is  fully 
captured  so  that  superposition  can  be  applied  is  to  cap  off  the  perforation  CAD  geometry  with  a  thin 
layer  of  material,  as  shown  in  figure  1 1 .  The  thickness  of  the  cap  layer  should  be  greater  than  the 
side  dimension  of  one  computational  element  so  that  the  geometry  can  be  resolved.  While  this 
approach  will  add  additional  material  to  the  grain  (and  thus,  alter  the  initial  volume  Vg),  by  allowing 
the  internal  feature  to  burn  immediately  at  the  beginning  of  the  analysis,  the  effects  of  the  initial 
excess  in  volume  become  negligible  as  the  burn  distance,  x,  exceeds  the  cap  thickness. 
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Figure  1 1 

Cross  section  of  single-perforation  grain  CAD  geometry  with  capped  perforation 

A  comparison  of  Vg  and  -dVg/dx  versus  x  for  a  single-perforation  grain,  with  numerical 
results  generated  with  and  without  the  superposition  technique,  to  the  exact  solution  is  presented  in 
figure  12  for  Nx  equals  301 .  A  cap  of  thickness  tcap  equals  0.0015  is  added  to  the  ends  of  the 
perforation  for  the  superposition  case,  which  is  slightly  larger  than  the  element  dimension 
Ax  equals  0.00134.  Slight  differences  in  Vg  and  dVg/dx  are  witnessed  at  the  onset  of  burning  due  to 
the  additional  cap  material;  however,  the  solutions  become  identical  at  x  equals  Ax/2,  when  the  cap 
is  completely  burnt  out.  As  x  increases,  the  two  solutions  differ  greatly  due  to  the  difference  in 
interface  treatment.  When  superposition  is  implemented,  numerical  results  conform  much  more 
accurately  to  the  analytical  predictions.  At  x  equals  0.12  (3.8%  of  initial  grain  remaining),  the  relative 
errors  are  23.9%  for  the  case  without  superposition  and  7.3%  for  the  case  with  superposition.  When 
superposition  is  implemented,  the  predicted  burnout  distance  decreases  from  0.151  to  0.127,  which 
is  within  2%  of  the  actual  value  of  0.125. 


Single  Perf  Grain  Regression  -  Superposition  Comparison  (Nx  =  301) 


- Volume,  No  Superposition - Volume,  Superposition  - Volume,  Analytical 

- dV/dx,  No  Superposition - dV/dx,  Superposition  - dV/dx,  Analytical 


6.00  E 


4.00  a 


> 

2.00  .E 


0.00  u 


Figure  12 

Grain  volume  and  change  in  volume  consumption  for  single-perforation  propellant  grain 

Seven-Perforation  Rod,  Further  Illustration  of  Superposition 

To  further  illustrate  the  importance  of  using  the  superposition  method,  consider  the  7- 
perforation  cylindrical  (rod)  propellant  grain  of  initial  diameter  Dg,  length  Lg,  perforation  diameter  Pg, 
and  web  thickness  wg,  as  shown  in  figure  13.  The  web  thickness  wg  of  a  7-perforation  rod  grain  (or 
the  distance  between  adjacent  perforations  or  between  perforations  and  outer  diameter)  is  defined 
as  wg  equals  Dg  -  3 Pg.  For  this  grain,  expressions  for  Vg  and  dVg/dx  for  a  given  burn  distance  x  < 
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2 

Wg/2  are  a  simple  extension  of  equation  12,  with  a  multiplier  of  seven  added  to  ( Pg  plus  2x)  in  the 
expression  for  Vg  to  account  for  the  additional  perforations.  For  x  >  wg/ 2,  however,  the  analytical 
expressions  become  more  complex,  as  the  intersections  of  the  perforations  both  with  the  outer 
boundary  and  each  other  cause  grain  slivering  (ref.  7),  and  thus,  are  not  repeated  here.  The 
superposition  method  discussed  in  the  previous  section  is  best  suited  for  predicting  burn  regression 
for  this  case,  as  it  will  be  able  to  better  capture  the  slivering  behavior  by  avoiding  diffused  interfaces 
and  small  geometries. 


Figure  13 

7-perforation  rod  grain  geometry 

Figure  14a  plots  a  comparison  of  Vg  and  -dVg/dx  profiles  versus  burn  distance  x  for  a  7- 
perforation  rod  grain  with  initial  Dg  equals  0.4,  Lg  equals  0.4,  Pg  equals  0.05,  and  wg  equals  0.0625 
generated  using  AGES  both  with  and  without  superposition  to  the  analytical  result.  The  numerical 
results  are  generated  on  a  cubic  domain  of  side  length  0.404  using  N  equals  201  elements  per  side. 
For  this  grain,  slivering  begins  at  x  equals  0.03125  and  is  indicated  by  the  sharp  decrease  in  dVg/dx 
after  the  web  has  burnt  out.  Total  burnout  occurs  at  x  equals  0.0512. 


(a) 


(b) 


Grain  volume  and  change  in  volume 
consumption 


Relative  percentage  error  in  grain  volume 
versus  burn  distance  for  7-perforation  rod 
propellant  grain 


Figure  14 

7-perforation  rod  superposition  comparisons 
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When  superposition  is  ignored,  significant  errors  in  the  numerical  results  occur  well  before 
slivering  begins,  and  continue  to  dramatically  increase  throughout  the  burn  duration.  As  seen  in 
figure  14b,  errors  of  over  15%  are  witnessed  at  the  onset  of  slivering,  and  the  total  burnout  distance 
over  predicts  the  analytical  by  27%.  When  superposition  is  included,  and  eight  separate  features 
(outer  boundary  and  7-perforations)  are  analyzed,  greater  agreement  between  numerical  and 
analytical  results  is  seen.  Errors  rise  to  ~7%  at  the  onset  of  slivering  but  remain  of  similar  magnitude 
until  the  grain  has  nearly  burnt  out  at  x  equals  0.0509,  which  is  within  1  %  of  the  analytical  value. 

Additive  Manufacturing,  “Fractal”  Grain 

As  mentioned  previously,  an  advantage  of  additive  manufacturing  is  the  ability  to  create 
complex  grain  geometries  containing  internal  structures  that  become  exposed  after  a  certain  portion 
of  the  initial  grain  burns  away.  Inclusion  of  complex  internal  structures  and  burn  surfaces  allow  for 
greater  control  and  fine-tuning  of  burn  rate.  A  prototypical  example  of  an  additively-manufactured 
propellant  grain  geometry  is  shown  in  figure  15,  which  is  referred  to  as  a  “fractal”  grain.  On  the 
exterior  (fig.  15a),  this  grain  appears  similar  to  a  cylindrical  rod  grain.  However,  inside  of  the  grain 
(fig.  15b  and  c),  a  large  internal  channel  that  spans  the  majority  of  the  grain  and  branches  out  in 
multiple  directions  exists  that  will  break  the  grain  apart  into  pie-shaped  fragments  and  dramatically 
increases  the  burning  surface  area.  Additionally,  as  shown  in  figure  15d,  24  smaller  internal 
subdivisions  exist  that  will  cause  grain  breakup  at  a  later  point  in  the  ballistic  cycle  to  further  increase 
grain  progressivity. 


(c)  (d) 

Angled  cross  section  showing  Angled  cross-section  showing 

larger  connected  internal  feature  smaller  internal  subdivisions 

Figure  15 

A  prototypical  example  of  an  additively-manufactured  propellant  grain  geometry 

In  the  previous  section,  non-physical  caps  were  placed  over  perforations  in  order  to  apply 
superposition  to  be  able  to  more  accurately  calculate  grain  regression  characteristics.  In  order  to 
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ensure  that  the  addition  of  cap  material  to  the  perforations  only  had  minimal  impact,  all  of  the 
internally  detected  features  were  allowed  to  burn  as  soon  as  the  numerical  analysis  began.  In  a 
case  such  as  this,  all  internal  features  should  begin  to  burn  when  first  intersected  by  another  burning 
interface,  whether  that  be  an  additional  internal  feature  or  regression  of  the  outer  grain  boundary. 
The  AGES  allows  for  a  user  input  to  dictate  whether  internal  features  are  allowed  to  burn 
immediately  or  based  on  intersection  with  another  feature  interface  boundary.  A  representative 
output  block  from  AGES  for  the  aforementioned  fractal  grain  is  shown  in  figure  16,  where  callouts  to 
the  burn  distance  at  which  specific  internal  features  begin  to  burn  is  explicitly  stated. 
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Figure  16 

Output  block  from  AGES  for  fractal  grain  with  N  equals  78  M  cells 

Figure  17  plots  grain  volume  and  change  in  volume  consumption  for  the  fractal  grain  shown 
in  figure  15,  generated  with  a  total  of  31 , 78,  and  128-million  computational  cells,  respectively.  For 
all  three  grid  resolutions,  the  numerically  calculated  initial  grain  volume  Vg  is  within  0.5%  of  the  actual 
value  with  slight  errors  occurring  due  to  curvature  approximation  using  cubic  elements.  As  the  grain 
begins  to  burn,  the  largest  internal  feature  starts  to  burn  at  an  analytical  burn  distance  of  x  equals 
0.19228,  which  is  the  shortest  linear  distance  between  the  surface  of  the  feature  and  the  outer 
boundary.  For  all  three  grid  resolutions,  the  numerically  calculated  burn  distances  of  x  equals  0.17 
to  0.18  are  accurate  to  within  the  edge  length  of  one  computational  cell,  which  is  the  maximum 
accuracy  that  can  be  achieved  using  the  current  method.  The  large  increase  in  total  burning  surface 
area  resulting  from  the  ignition  of  the  large  internal  feature  causes  a  spike  in  dVg/dx,  showcasing  the 
progressivity  of  the  grain. 
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Figure  17 

Grain  volume  and  change  in  volume  consumption  for  additively-manufactured  fractal  grain  with  total 

computational  cell  counts  of  31,  78,  and  128  million 

As  burning  continues,  a  smaller  spike  in  dVg/dx  is  witnessed  between  x  equals  3.05  to  3.15, 
resulting  from  the  ignition  of  the  24  smaller  internal  features  shown  in  figure  1 5d.  These  features 
begin  to  burn  when  intersected  by  the  already  burning  large  internal  feature.  The  initial  minimum 
distance  between  the  larger  internal  feature  and  the  24  smaller  features  is  x  equals  3,  and  thus, 
analytically,  the  burn  event  for  these  smaller  features  should  begin  at  x  equals  3.19228.  The 
numerically  predicted  burn  distances  where  the  internal  features  begin  to  burn  are  accurate  and 
within  1  to  5%. 

The  agreement  between  the  predicted  burn  distances  and  the  theoretical  values  obtained 
from  inspection  of  the  CAD  geometry  showcase  the  ability  of  AGES  to  provide  valuable  information 
about  the  burn  characteristics  of  arbitrarily-shaped,  complex  additive  manufactured  grains  without 
requiring  printing  and  testing  of  said  geometries. 


CONCLUSIONS 

The  ability  to  create  complex,  arbitrarily-shaped  propellant  grains  with  internal  features  using 
additive  manufacturing  processes  allows  for  greater  control  over  propellant  gas  generation  during 
munition  launch  cycle.  The  U.S.  Army  Armament  Research,  Development  and  Engineering  Center, 
Picatinny  Arsenal,  NJ,  grain  evaluation  software  (AGES)  was  developed  to  accurately  predict 
propellant  grain  regression  for  arbitrary  geometries  using  numerical  methods,  in  order  to  allow  for 
quick  analysis  of  grain  regression  characteristics  without  requiring  grain  printing  or  costly  burn 
testing.  The  AGES  has  been  validated  against  commonly-used  propellant  grain  geometries,  and  the 
capability  of  AGES  to  analyze  complex  geometries  with  intricate  internal  features  is  presented.  The 
AGES  is  government  owned,  developed,  and  maintained. 
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